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ABSTRACT 



The near proportionality between HI and dark matter in outer galactic disks prompted us to run N-body simulations of galactic disks 
in which the observed gas content is supplemented by a dark gas component representing between zero and five times the visible gas 
content. While adding baryons in the disk of galaxies may solve some issues, it poses the problem of disk stability. We show that 
the global stability is ensured if the ISM is multiphased, composed of two partially coupled phases, a visible warm gas phase and 
a weakly collisionless cold dark phase corresponding to a fraction of the unseen baryons. The phases are subject to stellar and UV 
background heating and gas cooling, and their transformation into each other is studied as a function of the coupling strength. This 
new model, which still possesses a dark matter halo, fits the rotation curves as well as the classical CDM halos, but is the only one to 
explain the existence of an open and contrasting spiral structure, as observed in the outer HI disks 
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1. Introduction 

The ACDM scenario encounters much success in reproducing 
the large scale structures of the Universe traced by th e Lyman-a 
forest and gravitational lensing (ISpringel et a DEqqI). However, 
at galactic scale, where the baryonic physics plays a major 
role, the ACDM scenario has sev eral well known problems . 
The dark matter cusp in galaxies (Blais-Ouellette etal] 12001 



the galax y formation process (Seiden et"anil984l) . However, as 
shown bv lvan den Boschl(l200lT the angular momentum conser- 
vation leads to the formation of too concentrated stellar disks, 
a real problem for the low surface brightness galaxies. On the 
other hand, the Hl-dark matter correlation may suggest that a 
large fraction of the dark matter lies in the disk of galaxies, fol- 
lowing the HI distribution. This physical link between HI and 
dark matter has been confi rmed by the baryonic Tu lly-Fischer 
relation of spiral galaxies dPfenniger & Revaz 20051) and mo re 
recently for extremely low mass galaxies ( Begum et al.ll2008l) . 

Dark matter in the disk of galaxies is now also sug- 
gested by dynamical arguments based on the asymmetries 
of galaxies, either in the plane of the disks or transverse 
to it. For example, the large spiral structure present in the 
very extended HI disk of NGC2 915 is supported by a 
quasi s elf-gravitating disk jBureau et a l. 1999; Masset & Bureaul 
l2003h . iRevaz & Pfennigej (12004') have shown that heavy disks 
are subject to vertical instabilities (also called bending in- 
stabilities) and may generate all types of observed warps; 
the common S-shaped but also the U-shaped and as ymmet- 
ric on es (Reshetnik ov & Co mbes 1999; Sanchez-Saav edra et al.l 
l2003h . Other arguments for the origin of warps have 
been proposed (iJiang & Binnevlll999l; IShen & SellwoodI 120061; 
[Weinberg & BUtzll2006l)T however none of them are able to ex- 
plain the three types of warps. 

Several works have tried to study the flattening of the 
Milky Way halo us ing potential tracers like dwarf galaxies 
dJohnston et alj2005h . From theses studies, the halo potential ap- 
pears to be nearly spherical, with an ellipticity of 0.9. However, 
dwarf galaxies trace the halo at distances larger than 20kpc 
where it is difficult to distinguish "classical" CDM disks with a 



de Blok & Bosma'2002';'S waters et al."2003';'Gentile et al.^2004 
2005; Spekkens et al. 2005; de Blok 2005; de Blok et al. 200i 
Spano et alj 120081) , and in particular in d warf irregulars, the 



high number of sinall systems orbiting halos dMoore et aljl999t 
iKlvpin et al.lll999l ; IStrigari et al.ll2007h . the low angular momen- 
tum problem at the origin of too small disks ( Nav arro & Benzl 
11991; Navarro & Steinmetz 1997; Kaufmann et al.l |2007|). as 
well as the difficulty of forming bulgeless disks ( iMaver et al.l 
120081) suggest that some physics is missing. 

A strong prediction of the ACDM scenario is that the so 
called "missing baryons" reside in a wa rm-hot gas phase in 
the over-dense cosmic filaments dCen & O striker 1999, 2006). 
However, there are now several theoretical and observational ar- 
guments that support the fact that galactic disks may be more 
massive than usually thought, containing a substantial fraction 
of these "missing baryons". 

It has been pointed out bv iBosmal (1 1 978i 1 1 98 lb tiiat in sam- 
ples of galaxies, the ratio between the dark matter and HI sur- 
face density is roughly co nstant well aft er the optical disk (se e 
also ICarignan etal]ll990l; lBroeilslll992l; iHoekstra et a"ni200l[) . 
This correlation may be a direct consequence of the conser- 
vation of the specific angular momentum of the gas during 
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slightly flattened halo from a heavy disk mo del, containing half 
of the halo mass in an extended flat disk (iRevaz & Pfennigen 
I2007h . 

A more accurate method is to trace the potential near the 
plan e of the disk, using f o r example the vertica l gas distribu- 
tion. iKalberlal ( l2003ll2004l) : lKalberla et al.l ( |2007|) have modeled 
the HI distribution out to a galactic radius of 40 kpc. Their self- 
consistent model is compatible with a self-gravitating dark mat- 
ter disk having a mass of 2 - 3 x 10" Mq. 

The presence of dark baryons in the disk of galaxies is 
reinforced by the numerous signs of recent star formation 
in the far outer disk of galaxies , coiTelated with HI gas in 
NGC 6946 dFerguson et al] Il998l) . N GC628 dLehevre & Rovl 



l205oh. M 1 01 ISmidi et al J (I2000I). M 3 1 I CuiUandre et al J (1200 1 



NGC 6822 dde Blok & Walteijl2003[) . M 83 dThilker et al. 



2005D. 



Correlations between young stars and HI far from the center re- 
veal that molecular gas, at the origin of a weak but existing star 
formation rate, must be present in abundance there, despite the 
lack of CO detection. 

Other arguments supporting the idea that dark matter could 
reside in the galactic disk in the for m of cold molecular hydr ogen 
(H2) have been widely discussed bv lPfenniger et alJdl994l) . 

Unfortunately, a direct detection of col d H2 in the outer disk 
of gal axies appears to be a very hard task dCombes & Pfennigeii 
I1997I) . As the H2 molecule is symmetrical, any electric dipole 
moment is canceled. The molecule may be detected in emission 
only by its quadrupole radiation, the lowest corresponding to a 
temperature of 512 K abo ve the fundamental state, mu ch above 
the 5 - 7 K expected from lPfenniger & CombesI dl994l) . A weak 
radiation of the H2 ultra-fine lines is expected at kilometer wave- 
lengths. However, since it is 10 orders of magnitude smaller than 
the HI line for the same density, its detection will not be possible 
in the near future. The absorption lines may be the best way to 
detect H2. But this method requires large statistics, as the filling 
factor of the gas is expected to be very low (< 1%). 

The best tracers of cold H2 may be the pure rotational lines 
of H2 (at 28, 17, 12 and 9 microns), which could be emitted 
by a few percent of the molecular gas, excited by intermittent 
turbulence (see for example Boulanger et al. 20081). 

Indirect detection by tracers may be prone to error. For ex- 
ample, CO traces the H2 but only for enriched gas and fails at 
large distances from galaxy centers. Moreover, it is impossi- 
ble to detect CO emission from a cloud at a temperature close 
to the background tem perature. The cold dust component de- 
tected by COBE/IRAS dReach et al]| 19951) is known to trace the 
cold H2. But as for the CO, it is limited to small galactic radii 
where the cold gas is still mixed with some dust. Spitzer mid- 
infrared observations have recently revealed that large quanti- 
ties of molecular hydrogen are not associate d with star forma- 
tion dAppleton et al.l 120061: [Ogle et al.ll2007l) . An unexpectedly 
large amount of 10^° M© of H2 is revealed o nly by its strong H 2 
emission lines in th e galaxy Zwicky 3146 dEgami et al.l 120061) . 
iGrenier et alJd2005l) showed that much gas in the solar neighbor- 
hood is revealed only by y-rays. Other indirect detections coul d 
be possible using micro-lensing events dPraind 1 9981; lFuxl2005h . 
However, this method has not been exploited up to now. 

As direct and indirect observational detection of cold gas is 
difficult, it is necessary to test the eff'ect a cold gas component 
would have on the global evolution of galactic disks. In this pa- 
per, we present new N-body simulations of galactic disks, where 
the observed gas content has been multiplied by a factor between 
3 and 5. In addition to this extra dark baryonic component, a 
non-baryonic, spheroidal pressure-supported dark halo contain- 
ing most of the large-scale dark mass is conserved. 



An important issue that ou r model aims to a nswer is the sta- 
bility question of heavy disks (lElmegreenll 1 997b . Our model as- 
sumes that the additional baryons h e in a very cold and clumpy 
phase dPfenniger & CombesI 1 19941) . partially dynamically de- 
coupled from the visible dissipative phase. We show that this 
phase can thus be less dissipative than the visible ISM and has 
larger velocity dispersions, so that the global disk stability is pre- 
served. A new numerical implementation of the cycling acting 
between these two phases is proposed. 

The secular evolution shows that the models with additional 
baryons are globally stable and share on average the same obser- 
vational properties as the "classical" COM disks. However, they 
give a natural explanation for the presence of contrasting spiral 
structures in the outskirts of HI disks which is diflicult to ex- 
plain when taking into account the self-gravity of the HI alone. 
In a forthcoming paper, we will show that this model also repro- 
duces th e puzzling dark matte r content present in debris from 
galaxies (iBournaud et al.ll2007l) . 

Our model is differen t from previous multiphase models 
(iSemehn & CombesI 12002"!: iHarfst et all 120061: iMerhn & Chiosil 
20071) in the sense that it does not compute a cycling between a 
cold-warm dissipative and hot medium, but between a very cold 
weakly collisional phase and the visible dissipative phase. 

The paper is organized as follows. Details of the multiphase 
model is given in Section |2] In Section [3] we briefly discuss the 
parameters used, and Section |4| is devoted to the galaxy model 
description. Section |5] compares the evolution of galaxy models 
with and without additional baryons and a short discussion. A 
summary is given in Section|6] 



2. The multiphase model 

2. 1 . The straightforward approach 

The straightforward approach when modelling the galactic ISM 
is to assume that gas behaves like an ideal, inviscid gas. The 
evolution of the specific energy u of the gas may be obtained by 
inserting the continuity equation into the first law of thermody- 



du P dp r(p) ■ 
dt p2 



Mp,t) 



(1) 



The first right hand side term corresponds to the adiabatic be- 
havior of the gas, while the second is responsible for the en- 
tropy variation. This latter reflects the non-adiabatic processes 
included through the heating and cooling function r(p) and 
A(p, T). At low temperature ISM (T < 10"^ K), the heating is 



dominated by the cosmic ray heating (iGoldsmith et al 



and by the ph oto-electric UV heating on grains ( WatsonI 



1969) 



19721: 



lDraindfT978h . As soon as the gas is enriched with metals, the 
coolin g function is dom inated by the CII, Sill, OI, Fell line emis- 
sions jMaio et al.l |200 7 t). Since the cooling function depends on 
the square of the gas density, the gas temperature is sensitive to 
it. 

However, the ISM is known to be strongly non homoge- 
neous down to very small scales, reaching densities higher than 
lO*" atoms X cm""* where the cooling time is very short, leading 
to equilibrium temperatures below 10 K. Such over-densities are 
unfortunately far from being resolved by galactic scale simula- 
tions and numerical simulations miss the associated low temper- 
atures. While being physically correct, Eq. [T]will then strongly 
bias the equilibrium temperature of the gas, because the esti- 
mated average density p poorly reflects the actual physics. 
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Fig. 1. Schematic representation of the ISM two levels system. 
The visible gas is labeled VG, the dark gas DG and the stars ST. 
The probability of transition between the visible and the dark gas 
depends on the UV flux generated by young stars (x-i,) and from 
extragalactic sources (x^xt) setting the characteristic transition 
times Tyg^dg and rag^vg- Stars are formed out of the visible gas 
only. 

2.2. Statistical approacli 

Instead of following the thermal specific energy of particles us- 
ing the biased Eq. [1] we propose a new statistical approach 
avoiding the problem of the density and temperature evaluation. 
The multiphase ISM is assumed to be a two level system (see 
Fig.lTJ with probabilities of transition depending only on the lo- 
cal excitation energy flux, that we call for short "UV flux", which 
is assumed to be the dominant heating process. 

The top level is populated by the well known observed 
dissipative gas detected by its CO, H2 or HI emission. This 
phase will be called the visible gas. The bottom level is popu- 
lated by undetected very cl umpy and cold gas as proposed by 
iPfenniger & CombesI (1 19941) . having temperature below the CO 
detection limit, at temperature equilibrium with the cosmic back- 
ground radiation. This gas results from the strong cooling that 
occurs in overdense regions. As this gas is missed by all tracers, 
we will call it the dark gas. 

According to the astrophysical literature, we have used here 
the word "gas". However, it is well known that the cold inter- 
stellar medium shows fractal properties which have been ob- 
served up to the instrument al capabilities, down to a few hun- 
dred AU ( lHeithausenll2"004l) . Such heterogeneous fluids clearly 
do not have the viscosity or other mean properties of smooth 
flows. In other fields, it would be called granular flows, for ex- 
ample. We are aware of the degree of simplification of our model 
compared to the complexity of the ISM. A more complex model 
should include the dark component, the CO-undetected metal- 
poor warmer H2 gas that may e xist in the outskirts of galactic 
disks dPapadopoulos et aLll2002l) . and possible effects related to 
phase transiti on and separation in the He-Ha mixture at very cold 
temperature (ISafa & Pfennigedl2008l) . 

The transition times between visible and dark gas is only 
dependent on the local UV flux which is the main heating 
process in the ISM. For the transition times, we have chosen the 
following simple relations: 



''"vg^dg 
T'dg^vg 



(2) 



where /3 and t are two free parameters, t gives the time scale of 
the transition, independently of the outer flux. f5 sets the flux de- 
pendency. For /3 -0, the transition does not depend on the outer 



flux. At equilibrium, these equations provide simple relations for 
the ratio between the visible gas density n^g, the dark gas density 
«dg and the total gas density «tg, as a function of the normalized 
local UV flux;if onlyQ: 



«vg 

«tg 
«dg 

«tg 
«dg 



2/3 



1 

-If) 



X 



(3) 



(4) 



(5) 



where in these latter equations, we have defined y as riya/ritg. 
The interpretation of these equations is the following. When the 
gas is subject to a strong UV field, the heating dominates and 
the probability of transition between the dark and the visible gas 
is high. The gas is then dominated by the visible phase. On the 
contrary, when the UV field is very low (as the case in the outer 
regions of the galactic disk), the cooling dominates and the prob- 
ability of a transition between the dark and the visible gas is low. 
The gas is then dominated by the dark phase. 

2.3. Visible and dark gas dynamics 

The visible phase is assumed to be dissipative, because it is 
more diffuse and collisional than the cold phase. Instead of us- 
ing the classical SPH approach, which makes the gas strongly 
collisional, we have preferred to use the sticky parti cle scheme 
(lBrahicll977l:ISchwardl98ll:ICombes & Gerinll985h which bet- 
ter simulates the clumpy ISM. In the following simulations, we 
have used /3i = 0.8 and ySt = 1 in order to strictly conserve angu- 
lar momentunfl The frequency of collisions between particles is 
set to be proportional to the local visible gas density. 

Also, as the very cold and clumpy dark gas does not radiate, 
it is expected to be weakly collisional. Its relaxation time being 
much longer than the dynamical time, we neglect here the effect 
of collisions. 



2.4. Star formation 

Stars that are assumed to be the main source of UV flux in the 
inner part of the galaxy may be formed out of the visible gas 
only. We have used a classical star formation recipe dKatz et al.l 
Il996h that reproduces well the Schmidt law. 



2.5. UVflux 

The normalized UV energy flux;^' is decomposed into two parts. 
A stellar radiation flux x* and an extragalactic background en- 
ergy flux;t'ext- 

2.5.1. Stellar UVflux 

The stellar normalized radiation flux is computed by summing 
the contribution of each star, assuming an appropriate (L/M) ra- 
tio. 



mi{L/M)i 



n + e- 



(6) 



' The solution is the one of a two level transition system with proba- 
bilities Pwc ~ l/T.rc and Pc„ ~ l/Tcf . 

- see iBournaud & CombesI ( l2002f) for details on the meaning of /S, 
and yS,. 
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j8 (L/M) (/ex, 

[erg/s/Me] [lO" Mo] 



1 1.58 X 10'- 1.8 X 10" 



[Gyr] 



10 



Table 1. Multiphase model parameters. 



where m, is the particle of mass i and r, its distance to the point 
where the flux is estimated, e is the gravitational softening, c the 
speed of light and C/g is the Habing UV energy density equal 
to 5.4 X 10 erg/cm^ (lHabing|ll968l) . The ratio (L/M) is set to 
1.58 X 10^' erg/s/Mg. This value is constant for all particles (in- 
dependently of the mass and type of star particles). The summa- 
tion in Eq.|6]is performed over all star particles and may be very 
time consuming in a large N-body system. We have computed it 
by taking advantage of its similarity to the summation of grav- 
itational forces which is computed using the treecode method. 
The time overhead for the computation of x* is flien negligible. 



2.5.2. Extragalactic UV flux 

The extragalactic UV background normalized radiation flux is 
assumed to be constant. 



(7) 



where the extragalactic UV background density is set to 2 x 
10"'^ erg/cm-^. 



3. Parameters 

Our multiphase model is based on 4 parameters. The external 
energy density flux f/ext, the UV light to mass ratio (L/M), the 
coefficient parameter and the time scale of transition t. In this 
paper, we have fixed the 4 parameters to constant values given in 
Table [1] The parameter is set to 1 for simplicity. The parame- 
ter (L/M) is chosen such that the normalized stellar UV flux^* 
is unity near 15 kpc when the stellar distribution at the origin of 
the UV flux corresponds to a realistic exponential disk, as will 
be presented in Section|4] According to Eq.|5] these values cor- 
respond to the radius where the dark component surface density 
is equivalent to the visible surface density. C/ext is chosen such 
that, in the absence of a stellar UV field (at a distance /? = oo 
from the galaxy center), the ratio «vg/«vg = Too where y^o has 
been set to 1/10. From Eq.[3]and setting = we get: 



t/ext = Ur 



1 



Tex 



(8) 



The time scale of the transition t is set to about 10 dynami- 
cal times so that the cycling between visible and dark gas is 
slow. Fig. |2] shows the total normalized flux x as well as the 
visible over total gas mass ratio y obtained with the parame- 
ters of Table [T] and the stellar distribution of Section |4] As an 
illustration, the flux is applied on a purely exponential bary- 
onic gaseous disk. The resulting surface densities for the visible 
and dark gaseous disk are displayed in the bottom of the figure. 
The exponential disk generates a strong UV flux at the center 
Consequently, the visible gas dominates there. On the contrary, 
in the outer part, the UV flux is dominated by the constant weak 
extragalactic flux, making the dark gas dominant. 




Radius [kpc] 



Fig. 2. Normalized UV flux x, visible to total gas density ratio 
■y and surface density of the gas, as a function of the galactic 
radius R. In the bottom graph, the red dashed line corresponds to 
the visible gas, the blue dotted line to the dark gas and the dark 
plain line to the total gas. 

4. Galaxy models 

We have used our multiphase model to simulate the evolution 
of Milky Way-like spiral galaxies. In order to understand the ef- 
fect of the multiphase model, we have compared two models, 
including additional dark baryons, with a reference "classical" 
model. In addition, we have build a fourth simulation including 
the perturbation of ACDM substructures. In the following sec- 
tions, we first present our reference model (also called N - I 
model), where no additional dark baryons have been added. In 
the next section (Section 14.1.2b we present how it is possible 
to add dark baryons in the reference model, in a way consistent 
with the rules proposed in Chapter|2]and by conserving the same 
observational properties of the reference model. The model in- 
cluding ACDM substructures will be presented in Section |4.1.3l 
Section 14.21 will discuss how initial velocities have been set in 
order to make the disk stable. 



4. 1 . Mass distribution 
4.1.1. The reference model 

Our reference model is designed to fit typical properties of giant 
Sbc galaxies with a flat extended rotation curve. It is initially 
composed of a bulge, an exponential stellar disk, a gas disk and 
a dark matter halo. 

1 . The bulge density profile is a simple Plummer model: 



<-Ij 



with a characteristic radius ry, - 1.4 kpc. 

The exponential stellar disk takes the usual form: 



p\R,z)oce 



-R/Hr ^-\z\IH, 



(9) 



(10) 



where the radial and vertical scale length are respectively 
//r = 4 and H, = 0.3 kpc. 

The dark matter is distributed in a Plummer model: 



'dm)' 



(11) 
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component 


model 


parameters 


mass 






[kpc] 


[10" Mo] 


bulge 


Plummer 


= 1.4 


0.28 


disk 


exponential 


//^ = 4,Hz = 0.3 


0.69 


halo 


Plummer 


'■dm = 30 


5.17 


gas 


Miyamoto-Nagai 


a = 30,b = 0.3 


0.21 



Table 2. Parameters for the reference giant Sbc model. 
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Fig. 3. Rotation curve of the reference model with the contri- 
bution of each component, the bulge, the exponential disk, the 
visible gas and the dark matter halo. 

with a characteristic radius ram = 30 kpc, ensuring a flat ro- 
tation curve up to 40 kpc. The model is truncated at 3 rdm 
(90 kpc). In order to avoid any perturbation of the halo on 
the disk due to imperfect equilibrium, we have set it as a 
rigid potential. We have checked that a live halo will not in- 
fluence our results. 
4. The choic e of the visible gas di stribution follows the obser- 
vations of lHoekstra et al.l (1200 lb . where the dark matter con- 
tribution to the rotation curve is a multiple of the contribu- 
tion of the gas. This is a chieved by distributing the gas in a 
Miyamoto-Nagai model dlVIivamoto & Nagalll 19751) : 



The proportionality between the Plummer halo (being in fact 
a subclass of a Miyamoto-Nagai model) and the disk is en- 
sured if the disk scale length is similar to the one of the halo 
: liR = '"dm (see Appendix [All. The vertical scale height is 
fixed to 0.3 kpc. We also have included a flaring of the visi- 
ble gas disk by multiplying the z coordinates of the particles 
by: 



exp(— ), 

where the transition radius Rf is set to 40 kpc. 



(13) 



Mass and parameters of the model are summarised in Table |2] 
and the corresponding rotation curve is displayed in Fig. [3] 

4.1 .2. Adding dark baryons in [he disk 

Here, we present a method that allows us to build new galaxy 
models containing additional dark baryons, but presenting sim- 



ilar observational properties as the reference model. The new 
models are build by following three rules: 

1 . The total rotation curve remains nearly unchanged with re- 
spect to the reference model. 

2. The surface density of the visible gas remains nearly un- 
changed with respect to the reference model. 

3. No dynamically significant dark matter is added in the cen- 
tral regions. 

In addition to these rules, the model, must also respect the ratio 
between the visible and additional dark gas that follows from the 
multiphase model proposed in Section|2] 

The third rule is derived frorti the luminosity MiUcy-Way 
model dBissantz & GerhardI 120021) cornbined with MACHO 
micro-lensing observations (see GerhardI ( l2006l) , and references 
therein for more details). We thus assume 1) no dark gas exists 
at the center of galaxies, and 2) the density of the total gas (tg) 
is equal to that of the visible gas (vg): 



p° = p^°, in the central regions. 



(14) 



This equation also ensures the second rule at the center 

In the outer part, the total gas mass distribution is constrained 
by the first rule. We can easily transfer mass from the halo to the 
disk without changing the rotation curve, using the equivalence 
in term of density of Eq. lA3] (via the Poisson equation). In that 
case, the density of the halo plus the total gas (htg) (including 
visible and dark gas) may be written as: 



hte 



(15) 



where the first left side term corresponds to the total gas density: 



tg . MN 



and the second left side term corresponds to the halo density: 

pk-(l-/«.o,..)- (17) 

In this equation, the mass of each component is given as an in- 
dex. Mhtg is the sum of the halo and gas mass of the reference 
model. Mtg = / Mhtg is the total gas mass and Mh = (1 - /) Mhtg 
is the halo mass. It is convenient to define / as: 



(12) f^N 



Mhtg ' 



(18) 



where gives the ratio between the total gas and the gas of the 
reference model. Its value is restricted to the range [0, Mhtg/Mg]. 
When = 1, Eq.[T6| and [T7| give the densities corresponding to 
the reference model. 

We can derive the total gas density all along the disk by com- 
bining Eq.[T4land[T6l 



p^R) = a{R)pl{R) + (1 - a{R)) pl,(R) 
= a{R) p^HR) + (I 



(19) 
(20) 



where we have introduced a new function a{R) taking values of 
1 at the center and in the outer part of the disk. 

We now introduce the multiphase model of Section |2| that 
places constraints on the ratio between visible and dark gas. 
Applying Eq. |4| to the total gas density ptg we can derive the 
density of the visible gas: 



p''S(R)^yiR)p\R) where 7(R) 



X' 



1 +x 



2/?' 



(21) 
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0.00 
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6.36 
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0.12 


0.42 


1.52 


6.27 


24 
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0.16 


0.66 


1.80 


6.12 


29 



Table 3. Parameters for the reference giant Sbc model. 



and where y(/?) (see Eq. |4|i follows from the stellar mass dis- 
tribution, the (L/M) ratio and the UV background. Combining 
Eq. |20] and |2T| allows us to determine the value of a{R) which 
completely defines the total gas distribution: 



1 



y(R) 



(22) 



In this equation, we force values of a to stay in the range [0, 1]. 
As we will see later, this restriction breaks the second rule in the 
outer part of the disk. Finally, the same flaring as for the visible 
gas (Eq.fTSll is applied to the total gas. 

Following this scheme, we have constructed two different 
models, with respectively N = 3 and N - 5. Applying the mul- 
tiphase model with parameters listed in Table [T] we can now 
determine the mass properties of the different components (vis- 
ible gas, dark gas, total baryons and dark halo) in these models, 
including the reference model (A^ = 1), given in Table[3] 

A larger value of increases the dark mass and inversely 
decreases the dark halo mass (the total mass remaining constant). 
The baryon fraction (Mb/Mtot) is then an increasing function of 
and grows from the universal baryonic fraction up to 30% 
for model N = 5. The total visible mass is not strictly identical 
between the three models. Its variation is mainly due to the the 
outer part where its surface density decreases with respect to the 
reference model. However, the visible gas lying below 20 kpc is 
similar for the 3 models (Mvg = 5.4 x lO^Mg). 

The gas and total baryon surface density is displayed in 
Fig. ID and |5] The outer regions (R > 20 kpc) of models N = 3 
and N - 5 are characterised by the presence of the dark gas 
which dominates the surface density of the baryons. In the far 
outer parts, the baryons are multiplied by N, with respect to the 
reference model. On the contrary, due to the drop in the dark 
gas surface density, the visible gas surface density of the inner 
regions is left unchanged with respect to the reference model 
(second rule). The presence of dark gas is only marginal inside 
10 kpc. 

Fig. |6] and |7] compare the contribution of all components to 
the rotation curve. In the models with additional baryons, the de- 
creasing contribution of the dark halo mass in the outer part is 
compensated by the dark gas, ensuring a flat rotation curve up to 
50 kpc. The bump in this component, appearing below 25 kpc, is 
due to the fact that, as no dark gas resides in the central regions, 
the square velocity of this component alone is negative. As in 
Fig. |6]and|7] we have plotted the absolute value of the velocity, 
the imaginary part appears as positive. Except around the tran- 
sition radius of 20 kpc, the total rotation curve of models N - 3 
and N = 5 remains nearly unchanged with respect to the refer- 
ence model (first rule). 

4.1.3. ACDM substructures 

In addition to the three previous models, we have built a fourth 
model called the = 1 + s model based on the N - I models 
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Fig. 4. Comparison of the surface density of baryons between 
model N = I and model N - 3. The total baryons (including 
stars) are traced in black. The red line represents the visible gas 
while the blue line (only present for the N = 3 model) falling 
towards the center corresponds to the dark gas. 
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Fig. 5. Same figure as Fig. |4]but for model N - 5. 



but taking into account a sample of dark matter satellites pre- 
dicted by the ACDM cosmology, orbiting around the disk and 
interacting with it. 

The purpose of this model is to study the effect of ACDM 
satellites on the spiral structure of the disk, and to compare it 
with the effect of the additional baryons taken into account in 
models A^ = 3 and N - 5. 

The ACDM sat ellites have been adde d following the tech- 
nique described by Gauthier et al.l (l2006h based on the distri- 
bution function of iGao et all ( |2004|) . in agreement with recent 
ACDM simulations. 100 satellites extending up to 250 kpc, with 
masses between 1.5 x lO""* and 0.02 mass of the galaxy of the 
reference model (6.36 x 10" Mq), have been ad ded, represent- 
ing 10 % of the total galactic mass. Contrary to iGauthier et al.l 
(|2006|) . in our simulations, the substructures are Plummer 
masses with a softening fixed to 0.25 kpc. Fig.[8]shows the sur- 
face density of model A^ = 1 -i- s including the ACDM satellites. 
The darkness of the satellites scales as a function of the mass. 
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Fig. 6. Comparison of the rotation curves of model N = I and 
model N = 3 and the contribution of each component.) 
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Fig. 7. Same figure as Fig.|6]but for model N = 5. 




Fig. 8. Model = 1 + s : 100 ACDM satellites around the 
gaseous disk of the N - I model. The images have been take 
after t = 2.2 Gyr of evolution. The box size is 400 x 400 kpc. 



4.2. Velocity distribution 

The velocity distributio n s are computed following the method 
proposed by lHernquis"3 ( Il993l) . For the spherical distribution, 
the bulge, the velocity dispersion is assumed to be isotropic 



and is derived from the se cond moment of the Jeans equation 
dBinney & Tremain3ll987h which is then: 



1 r°° 

Phir) Jr 



(23) 



For the axisymmetric components (stellar disk and gas 
disk), we first compute the vertical velocity dispersion cr-, 
by satisfying the Jeans equation in cylindrical coordinates 
(iBinney & Tremain3ll987h separately for each component /: 



1 r'" 

— dz'piiz 

^iiz) Jz 



(24) 



For the stellar disk and for the bulge, the radial velocity disper- 
sion (Tr is fixed by imposing the value of the Safronov-Toomre 
parameter Qto I, ensuring the horizontal stability of the disk. Q 
is defined by: 



Q 



KCTr 



3.36 G (Erf + 2;„) 



(25) 



where k is the radial epicyclic frequency. For the gas, according 
the the HI observations, ctr is fixed to a value of lOkm/s. The 
tangent ial velocity cr^ is derive d from the epicycle approxima- 
tion (BinniyI&^maini[T987l p. 125): 



a] ,2 R-'dR[R'dR^) 



z=0 



o-l 4Q2 AR-^ (5«(D),,=o 



(26) 



where the rotation frequency Q is determined from the poten- 
tial /^-derivatives as indicated in Eq. ( |26] |. Finally, the mean 
azimuthal velocity is deduced from the second moment of the 
Jeans equation: 

<y2> = R dR(l>(R, z) + o-l-o-l + ^Or {picrl) , (27) 
where we have assumed that the term {vrv^} - 0. 

4.2.1 . Velocity distribution for tine visible and dark gas 

The radial velocity dispersions of the visible gas is set to be con- 
stant at lOkm/s as observed: 



o-R^„ - lOkm/s. 



(28) 



On the contrary, in order to ensure stability, the radial veloc- 
ity dispersions of the dark gas is set by imposing the Savronov- 
Toomre parameter to be 1 : 



0"Rdg 



3.36 G (Erf + E„)e 



with Q- I. 



(29) 



The velocity dispersions and the mean azimuthal velocity are 
then used to distribute the model particles in the velocity space. 

5. Model evolution 

The multiphase mo del has been im plemented on the parallel 
Tree code Gadget2 (ISpringelll2005h . The models = 1, = 
1 -I- s, = 3 and = 5 contain respectively 319315, 319415, 
535769 and 731129 particles. The mass of the gas particles is 
constant between the three models. The softening length is set 
to 250 pc. All simulations have been run between and 4.5 Gyr. 
For simplicity, the feedback from supernovae has been turned 
off^. 
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Fig. 9. Comparison of the rotation curve of the 3 models at r = 
4.5 Gyr. 
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Fig. 11. Comparison of the surface density of different compo- 
nents at f = 4.5 Gyr. 
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Fig. 10. Surface density map of the stellar disk of the 3 models 
(A? = 1, A? = 3 and A? = 5, from left to right) at t = 4.5 Gyr. The 
box size of each image is 60 x 60 kpc. 



5.1. Global properties 

Despite their different dark matter content, after 4.5 Gyr, the 
three models {N = I, N = 3 and N - 5) still share similar 
global properties. 

In Fig. |9] we compare the total velocity curve of the three 
models, after 4.5 Gyr of evolution. In the central part, they all 
show a bump corresponding to the presence of a bar In the 
outer part, all curves converge to the same values. The main 
differences occur around the transition radius (~ 10 - 20 kpc), 
where models with additional baryons have slightly lower val- 
ues. However, these differences are simply a relic of the differ- 
ences existing in the initial conditions (see Fig. |6] and |7|i. The 
similarity of the rotation curves means that the three models 
share the same radial potential dependency, having similar hor- 
izontal epicyclic frequencies {k and Q.). The pattern speed, the 
shape and the extension of the bar is also identical (see Fig.fTOli. 

The surface density of the three models is compared in 
Fig. [TT] For radius smaller than 30 kpc, the curves corre- 
sponding to the visible gas surface density are superimposed. 
Differentiating these three curves is only possible in the far outer 
part, as was the case for the initial conditions (see Fig. |4] and 
|5]l. It is thus important to notice that the rotation curve and the 
azimuthal averaged surface density are not sufficient to easily 
distinguish between these three models, even after 4.5 Gyr of 
evolution. 

In FiglT2]we compare the scale height of the disk of the four 
models at f = 0.9 and at f = 4.1 Gyr. The continuous lines corre- 
spond to the visible gas while the dotted ones to the dark gas. At 
t - 0.9 Gyr, no difference exists between the four models. After 




20 30 
Radius [kpc] 

Fig. 12. Disk scale height as a function of the radius, for the four 
models at f = 0.9 and f = 4. 1 Gyr. The lines corresponds to the 
visible gas while the dotted lines correspond to the dark gas. 



4 Gyr of evolution the visible gas in the = 1h-s, A^ = 3 and 
N - 5 models has been slightly heated and presents higher scale 
heights. In the case of model = 1 -i- 5, the increase of the scale 
height is attributed to the perturbation of the ACDM satellites 
(Fo nt et a l. 2001; Dubinski et al. 2008), while for models A^ = 3 
and A^ = 5 it is due to the coupling with the collisionless dark 
gas that has a higher scale height (dotted lines). From an obser- 
vational point of view, models A^ = 1 + A^ = 3 and N = 5 will 
be very similar. 

5.2. Disk stability and spiral structures 

The main differences between the models with additional dark 
baryons and the "classical model" appear when comparing the 
spiral structure of the visible gas. The surface density evolution 
of the three models plus the one with the ACDM satellites is 
displayed in Fig. [T3lfT6] In the two latter plots, in addition to the 
visible gas surface density, the dark gas surface density is also 
displayed. These plots emphasize two important points. First, 
despite the presence of larger amount of baryons in the outer 
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Fig. 15. Evolution of model N - 3 between and 4.5 Gyr. The first and third columns represent the visible gas surface density while 
the second and fourth represent the dark gas surface density. The box size is 200 x 200 kpc. 
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Fig. 16. Evolution of model N - 5 between and 4.5 Gyr. The first and third columns represent the visible gas surface density while 
the second and fourth represent the dark gas surface density. The box size is 200 x 200 kpc. 
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parts, the disks of models N - 3 and = 5 are globally sta- 
ble over the 4.5 simulated Gyrs. Secondly, spiral structures have 
very different patterns and may be used to distinguish between 
models with and without additional baryons. 

At all times during the evolution, the surface density of the 
dark gas decreases at the center, preserving a hole. This hole 
simply means that no additional dark matter is expected at the 
center, even after 4.5 Gyr of evolution. 

5.2.1. Disk stability 

The stability of models N - 3 and = 5 is ensured by the 
presence of the two partially dynamically decoupled phases. The 
dark phase has higher velocity dispersions than the visible gas. 
The radial velocity dispersions of model N = 3 are displayed in 
Fig. [17] Because the dark gas is quasi non-collisional, its radial 
velocity dispersion is much larger than the one of the visible gas 
of the reference model (values dropping to lOkm/s in the outer 
part.). These higher values ensure the stability of the total bary- 
onic disk. However, our multiphase model fails to reproduce the 
low velocity dispersions expected in the visible gas. While be- 
ing clearly decoupled from the dark component, the gas velocity 
dispersion of model = 3 is higher than the one of the reference 
model, especially in the outer part. This is the result of the weak 
coupling due to the cycling between the dark and the visible gas. 
The model could be improved in the future, by simply increasing 
the stickiness of particles in the outer part. 

As discussed in AppendixiBl if one assumes that the dark gas 
is as dissipative as the visible gas, the resulting velocity disper- 
sions of the total gas is no longer large enough to ensure the disk 
stability. In that case, the disk breaks and forms small clumps 
(Fig.lO). 

Since our disk model is made of three dynamically de- 
coupled components (the stellar disk, the visible disk and the 
dark disk) having distinct velocity dispersions, a precise study 
of its stability would require to use a multi-component stabil- 
ity criterion. The sta bility of multi-compone nt galactic disks 
has been discussed bv lJog & Solomoril ( Il984bllal) and the equiv- 
alent o£the Savronov-Toomre parameter may be computed 
(lJoglll996h . While this criterion is theoretically valid for an n- 
components system the computational procedure described by 
Jod ( 19961) is only v alid f or a two-componen ts system (See also 
Wang & Silkl (11994 and lElmegreenI [l995l) ). We have instead 
used a Savronov-Toomre (Q\-c) assuming a single component. 
Formally, this 2i-c is defined by: 



3.36G£disk 



(30) 



where Zjisk and cr^^disk are computed over the visible gas, dark 
gas and stars particles, while k results from the total potential (in- 
cluding halo dark matter and disk dark baryons). The left panels 
of Fig.[T8]shows the parameter gi-c for the four models at differ- 
ent times. It is clear that Qi-c is always larger than 1, explaining 
the stability observed in Fig.[T5]and[T6l 

5.3. Spiral structure and swing amplification 

While being larger than 1, the two models with additional 
baryons are characterized by a nearly constant Qi-c from 10 to 
50 kpc, with values between 2 and 3 for model N = 3 and around 
2 for model N - 5. With these rather low values, the disks are 
only quasi-stable. Contrasting and open spiral structures are con- 
tinuously generated up to the end of the dark gas disk at 100 kpc. 
These spirals are thus naturally present in the visible gas up 
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Fig. 13. Evolution of model N = I between and 4.5 Gyr. Each 
image represents the surface density of the visible gas. The box 
size is 200 x 200 kpc. 
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Fig. 14. Evolution of model = 1 including 100 ACDM sub- 
structures. Each image represents the surface density of the vis- 
ible gas. The box size is 200 x 200 kpc. 
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Fig. 17. Radial velocity dispersion at f = 4.5 Gyr of the visi- 
ble (red) and dark (blue) gas components of model N - 3. The 
dotted line corresponds to the radial velocity dispersion of the 
visible gas of mode N - I. 



to its end, where its surface density drops, in agreement with 
most HI ob served disks. See for e xample the impressive case of 
NGC 6946 (iBoomsma et al.ll2008l) . 

On the contrary, model N - I has a higher Qi-c, increas- 
ing well above 3 for R > 30 kpc. Consequently, the outer disk is 
more stable, preventing the formation of spirals, as observed in 
Fig. [T3] This difference is the key point that allow us to distin- 
guish between models with and without additional baryons. 

In order to improve our understanding of the spiral structure, 
we have computed polar maps {{R - 0)-plot) of the visible gas 
surface density (from 7? = to = 50 kpc) at different times 
(Fig.[T9l). These maps can be, fo r example, compa r ed to t he HI 
observation in M83 (Fig. 12 of ICrosthwaite et al.l ( |2002|) ). We 
then performed a Fourier decomposition of those maps, for each 
radius. The spiral structure can be represented by each dominant 
azimuthal modes for a given radius, i.e. the mode with the largest 
amplitude, excluding the m - mode. Fig. |20] displays the am- 
plitude of the dominant azimuthal modes found for each plot of 
Fig. [19] For a direct comparison between the models, in all plots, 
the signal amplitude is coded with the same colors. In polar 
plots, a regular spiral structure of mode m appears as 2 m inclined 
parallel linefl where the inclination (that may vary with radius) 
gives the pitch angle ; of the arm (; = 7r/2 - arctan(p), where 
p is the slope of the lines). Such regular features are present for 
the N - 3 and N - 5 models, indicating the presence of rather 
open spiral arms (large that may be followed up to 50 kpc at 
nearly all times. The number of arms is not constant and oscil- 
lates between 3 and 8 (model N = 3) and 3 and 6 (model - 5), 
usually increasing with increasing radius (lBottemall2003h . Such 
large features are not observed for the N - 1 model. The spiral 
arms disappear around R -25 kpc. 

Spiral arms in a di sk may be the result of 

'swing ampUficatio n 'GGoldreich & Lvnden-Belll 119651; 
lJulian & Toomrd 119661; iToomrd 11981 ). a mechanism that 



locally enhances the self-gravity and leads to the amplification 
of a small perturbation in a differentially rotating structure. The 
amplification results from the resonance between the epicyclic 
motions of stars (or gas) and the rate of change of the pitch 



^ Here, a line is either a maximum or a minimum. 
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Fig. 18. Stability of models = 1, = 3 and = 5 at dif- 
ferent times (0, 0.9, 1.8, 2.3, 3.2 and 4.1 Gyr). The left column 
show the Toomre parameter (Eq. [30l l. The right column shows 
the mode corresponding to the swing amplification parameter 
X = 3. Modes higher than the latter may be amplified while 
lower modes are stable. In both columns, the dotted lines cor- 
respond to the values computed without taking into account the 
dark gas. The corresponding color code for all times are given in 
the upper right panel. 



angle of a density wave, during its conversion from leading to 
trailing ( iBinnev & Tremaing Il987h . The swing amplification 
mechanism may be quantified by the parameter X, the ratio of 
the perturbation wavelength and a critical wavelength: 



X = A/A^ 



where 



(31) 



Since the wavelength A of a spiral pattern of mode m in a disk is 
2nR/m, X becomes; 



XiR) 



(32) 



Similarly to the Qi-c computation, we have computed the swing 
amplification parameter X, taking into account only the surface 
density of the disk (E = Sdisk)- As it is too dynamically hot, at 
first order, the dark matter is not sensitive to the density pertur- 
bation and its effect is not included in the computation of X. If 
we assume tha t the swing arnplifica tion is efficient for X < 3 
(lToomrelll98lt lAthanassoulal 11984 Fig. 7) 0, we can reverse 
Eq.[32]and compute the value of the critical mode m(X - 3); 



m{X = 3) = 



Rk^ 



(33) 



disk 



In fact the exact critical value may vary, depending on the parameter 
and on the exact velocity curve through the shear parameter T = 
~n Jr ^ which in our case is nearly equal to 1. 










Fig. 19. Evolution of the spiral structure in a (/? - 0)-plot of the 
four models. Each image represents the surface density of the 
visible gas. The abcissa corresponds to the radius R (from to 
50 kpc) while the ordinate to the azimuthal angle (from to 27r). 
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Fig. 20. Evolution of the spiral structure in a (/? - 0)-plot of 
the four models. Each image represent the dominant mode of 



The meaning of the critical mode is the following: modes higher 
than m{X - 3) may be amplified by the swing amplification 
mechanism. On the contrary, the ones above m{X = 3) are stable 
and thus no spiral pattern with m < m{X - 3) is observed. 

The right panels of Fig. [T8]displays the critical modes at dif- 
ferent times. For the N = I model, m{X - 3) increases very 
quickly, being above 6 for R > 20 kpc. At larger radii, only large 
m modes, corresponding to noise, may be amplified by the self- 
gravity of the disk. No spiral pattern is expected, in agreement 
with Fig. |20]and[l3] 

When increasing Sjisk by including additional baryons, self- 
gravity in increased and lower modes may be amplified at larger 
radius. For example, m = 6 modes may still be amplified up 
to /? = 35 kpc in model N = 3 and up to /? = 50 kpc in model 
N = 5. This explains the patterns present in Fig.|20]and the spiral 
structure of Fig. [15] and Fig. [T6] 

The dotted lines of Fig.[T8](both in the left and right panels) 
correspond to the Qi-c and m{X - 3) values as would be de- 
duced by an observer ignoring the presence of the dark gas and 
taking into account only the visible gas and stars in Eg. l30land 
[33] Clearly, gi-c is well above 4 for a radius larger than 20 kpc, 
indicating that the disk should be stable. The m(X = 3) param- 
eter implies that modes with low perturbation should be prohib- 
ited by the swing amplification theory, in contradiction to the ob- 
served spiral structure. This point may explain why spiral struc- 
tures are present in galaxies where the Q\-c (or X) parameter 
deduced from the visible component is well above 5 (resp. 3), as 
is the case, for example, for the galaxy NGC 2915 (iBureau et al.l 
il999h . 

5.3.1 . The effect of the ACDM satellites 

We can now ask if the perturbations generated by the ACDM 
satellites are able to reproduce the large scale spiral patterns of 
galaxies, comparable with the ones obtained when including ad- 
ditional baryons. 

The effect of the satellites may be seen by comparing the 
surface density of the visible gas disk of model N - I (Fig.[T3]l 
and model N = 1 + s (Fig. [T4l i. Perturbations in the far outer 
part of the disk are visible between 2.3 and 4.1 Gyr and appear 
as asymmetric and winding spiral arms. These perturbations are 
better traced in Fig.|20l where the same Fourier analysis has been 
performed for model N = I + s. There, the effect of the satellite 
is already perceptible at r = 1 .4 Gyr and R ~ 30 kpc. Between 
f = 1.4 and t = 2.3 Gyr an m - 2 mode develops and disappears 
afterwards. The large slope of the modes indicates a very small 
pitch angle. Another m - 2 mode appears briefly between t - 
2.3 and t - 2.7 Gyr (R > 40 kpc). In this plot, no clear dominant 
mode is seen afterwards, between 25 and 50 kpc. However, as 
the Fourier decomposition was performed within R - 50 kpc, 
we are missing the spirals observed at larger radii in Fig. [14] at 
t = 3.6 kpc. 

The comparison between the spirals resulting from the self- 
gravity of the disk or the ones formed by satellites perturbations 
reveals two important differences: 

1. As it does not result from a gravity wave, a spiral arm 
generated by the perturbation due to a satellite is affected 
by the differential rotation of the disk. Consequently the 
arm winds quickly (the pitch angle is small) and disap- 
pears in a d ynamical time (winding problem, see for example 
iBinnev & T remaine ( 19^%). 

2. Spiral patterns generated by the swing amplification are 
globally continuous along the disk. We have seen for ex- 
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ample that the dominant modes slowly decrease outwards. 
As they result from a local perturbation, asymmetries and 
strong discontinuities exist in the pattern of spirals induced 
by satellites. An example is seen in Fig.|20]at f = 2.7 Gyr. 
No dominant mode exists between 25 and 40 kpc while an 
m = 2 mode is present at larger radius, disconnected from 
the rest of the disk. 

In conclusion, the ACDM satellites are able to perturb a disk 
and induce spiral arms. However, these arms wind up and disa- 
pear quickly. Even if they are periodically induced by satellites, 
they are unable to reproduce a regular and extended spiral pat- 
terns. 

According to the theory of spiral structure formation, the fact 
that ACDM induces only localised and decorrelated perturba- 
tions is not surprising. Indeed, if satellites act as triggers for 
disk instabilities, the subsequent amplification of the instabili- 
ties leading to the formation of a large scale pattern is deter- 
mined by the property of the disk itself: its self-gravity (surface 
density), its velocity curve and its velocity dispersions, the three 
quantities involved in the Qi-c and X computation. Without the 
interplay of these essential ingredients, only discontinuous local 
structures, resulting from the random impact of the satellites on 
the disk, will be observed. 



6. Discussion and conclusion 

Using a new self-consistent N-body multiphase model, we 
have shown that the hypothesis where galaxies are assumed 
to have an additional dark baryonic component is in agree- 
ment with spiral galaxy properties. The key element of the 
model is to assume that the ISM is composed out of two 
partially dynamically decoupled phases, the observed dissi- 
pative gas phase and a dark, very cold, clumpy and weakly 
coUisional phase. Motivations for invoking a similar com- 
ponent made of cold invisible gas, like the disk-halo con- 
spiracy, the Hl-dark matter proportionality, extreme scatter- 
ing and microlensing events, the hydrostatic equilibrium of 
the Galaxy or the survival of small molecular structures with- 
out s hielding, have been given many times in the literature 
(e.g.. iPfenniger & Combesril994t iHenriksen & Widrowl 119951: 
Gerhard & Silk" 1996V Walker & Wardle"l998\l Kalberla & KerpI 
1998 : Kerins et al. 2002; Heithausen 2002, 200^. 

An original scheme, designed to overcome numerical lim- 
itations, has been proposed to compute the cycling between 
these two phases when subject to heating and cooling. From 
this scheme, we have shown that realistic self-consistent N-body 
models of spiral galaxies containing an additional baryonic dark 
matter content may be constructed. These models share similar 
observational properties with classical CDM disks, like the rota- 
tion curve and visible gas surface density and in that sense, are 
observationally similar. 

The main result of this work is that, despite having more 
mass in the disk, these systems are globally stable, the stability 
being ensured by the larger velocity dispersion of the dark gas 
that dominates the gravity in the outer part of the disk. In addi- 
tion, the enhanced self-gravity of these disks, due to the presence 
of the dark clumpy gas, makes them more prone to form spirals 
extending up to 100 kpc in the dark gas. The spiral structure is 
revealed by the HI that acts as a tracer of the dark gas, up to the 
radius where its surface density becomes so small that the hydro- 
gen is fully ionized and hardly observable. This gives a natural 
solution of the numerous observations of HI unstable disks that 



are difficult to e xplain by the self-g r avity of the HI disk alon e 
(see for example lBureau et alJ ( ll999l) : lMasset & Bureaul ( l2003l) ). 

Depending on the theory of spiral structure (the local sta- 
bility of differentially rotating disks and the swing amplifica- 
tion mechanism) our results depend on the three quantities in- 
volved in the Toomre parameters: the disk self-gravity by its sur- 
face density, its rotation curve and its radial velocity dispersion. 
Consequently, these results will not be affected by a different 
parametrisation of the stellar disk or the dark halo, as long as a 
similar rotation curve is reproduced. 

We have seen that in our models the velocity dispersion 
of the collisionless or weakly collisional components increases 
with time, similarly to what stars are known to do in the 
Milky Way in the so-called Wielen's diffusion. This is fully ex- 
pected from our general understanding of the evolution of self- 
gravitating disks of collisionless particles, where the effects of 
spiral arms and especially a bar are sufficient to cause such a ra- 
di al heating. A more detailed discussion and references are given 
in lSellwood & Binnevl(l2002l) . 

While the cold gas phase is probably weakly collisional, we 
have treated it as strictly collisionless. How will the spiral mor- 
phology then be modified by introducing a weak dissipation in 
this component ? In Appendix IB] we present the case where the 
dark component is as dissipative as the observed gas. In this ex- 
treme case, the disk is strongly unstable and fragments as a con- 
sequence of the Jeans instability. At some point, decreasing the 
dissipation of the dark gas strengthens the global stability, while 
reinforcing the spiral arm contrast. The dynamical coupling of 
the two phases through the cycling is as important as the dissi- 
pation. A very short transition timescale t leads to a strong cou- 
pling. Additional simulations show that in that case, the velocity 
dispersions of the two components are nearly similar, the one 
of the observed gas being higher. Consequently, the spiral struc- 
tures are less contrasting. On the contrary, if the two phases are 
decoupled, the observed gas velocity dispersion is smaller, and 
the arm contrast higher. In both cases however, the large scale 
spiral pattern remains similar. 

Additional simulations have investigated the effect of 
ACDM substructures on the outer HI disk. Satellites only gen- 
erate local, winding and short-lived arms. Because the HI disk 
is not sufficiently self-gravitating, its response is too small to 
a mplify the perturbation s. In that sense, we confirm the results 
of iDubinski etai] ( |2008|) . The tidal effects of the satellites are 
generally small and so are not responsible for large scale spiral 
patterns. 

If models of galactic disks with additional dark baryons are 
in agreement with observations, we must ask whether they are 
in agre ement with the A CDM scenario. From the cosmic baryon 
budget dFukugita & Peebles 2004) . we know that nearly 60% of 
the baryons predicted from primordial nucleosynthesis are not 
observed. Multiplying the galactic baryon content (about 8%) by 
a factor of 2, as it is assumed in our most massive N - 5 model, 
will still be in agreement with the baryon budget. From a galactic 
point of view, multiplying the baryon content by a factor 2 is well 
inside the observed uncertainties. For a circular velocity curve of 
about 250km/s typical of our models, the baryon content varies 
by more than a factor of 5 (see Fig. 1 of Maver & Moore ( 2004)). 
Assuming a corresponding virial mass of 10'^ Mo the cosmic 
baryon fraction may vary around 10 and 50%. With a baryonic 
fraction of 29%, our N - 5 model corresponds to a large but 
realistic value. 

Numerical simulations of the formation of large scale struc- 
tures predict that the "missing baryons" reside in a warm-hot 
gas phase in the over-dense cosmic filaments dCen & Ostrikeii 
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Il999ll2006l) . However, there are now also signs of accretion of 
cold gas during the build up of galactic disks (iKeres et"al1l2005l) . 
As discussed in Section [2] a low density resolution leads to an 
under-estimate of the gas cooling, missing over-dense regions 
present in an inhomogeneous ISM, where the cooling time of 
the gas is very short. Since the assumed amount of dark baryons 
is not in contradiction with the cosmic baryon budget, we need 
to wait for future works to see if the correct treatment of inho- 
mogeneous ISM, including shocks, is able to cool 10% of the 
"missing baryons" during the hierarchical structure formation, 
making the present model consistent with the ACDM scenario. 

In a companion paper, we have shown that our model also 
explains the puzzling presence of dark mat ter in the collisional 
debris from galaxies jBournaud et al.ll2b07h . 

The remaining questions are obviously why the cold dark 
phase is invisible and what the physics ruling it might consist of. 
Since star formation is also a poorly understood process, but we 
know that stars do form from the coldest observable molecular 
gas, the question of why stars are not supposed to form in the 
coldest form of our invoked dark gas remains open. Suggestions 
of why s tar formation may not always start in cold gas have been 
given in iPfenniger & CombesI (Il994b . Essentially, if molecular 
gas fragmentation goes down to sub-stellar mass clumps stay- 
ing cold, these clumps are unable to free nuclear energy and 
become stars. Since they are close to the 3 K background, they 
almost do not radiate. Because of their negative specific capac- 
ity, such self-gravitating clumps do not increase their tempera- 
ture when subject to external heating, but stay cold while evap- 
orating, a sometimes overlooked feature of self-gravity. Many 
earlier works have discussed these issues, for example how 
cold molecular hydrogem may clu mp in dense structures of so- 
lar system size but stay undetected (ICombes & Pfennigeiill997t 
IPfenniger & CombesI fl994t iPfenniged |2004 . However, much 
more work remains to understand the precise physics of very 
cold gas (T < 10 K) in low excitation regions, like the outer 
galactic disks, but also in planetary nebulae and star forming re- 
gions where observed dense cold gas clumps self-shield from 
outer radiation and reach sub-lOK temperatures. It is gener- 
ally known that H2 condenses in solid form even at interstellar 
pressures at temperatures close to the 3 K cosmic background. 
A phase transition means that a richer physics must be ex- 
pected. But since H2 is usually mixed with about 10% of he- 
lium in mass, it was for many years not clear how to describe the 
equation of state of this mixture in such conditions. Recently, 
ISafa & PfeniiiieJ (|2008|) succeeded in describing this mixture 
for astrophysically interesting conditions with chemo-physical 
methods, reproducing its main characteristics like the critical 
point and the condensation curve, as well as predicting the con- 
ditions of He-H2 separation. Work is in progress to apply these 
results to the cold interstellar gas. 
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Appendix A: The Miyamoto-Nagai potential : 
another parametrization 



The Miyamoto-Nagai potential (fM ivamoto & Nagail 1 1975!) in 
cylindrical coordinates is usually written as: 



GM 



+ (fl + V?TF)^ 



(A.l) 



where M is the total mass and a and b are two parameters. It is 
useful to parametrize the Miyamoto-Nagai potential, introducing 
the horizontal scale length Iir and the vertical scale length h~ 
defined by: 



liR- a + b and h, - b. 
The potential then becomes: 
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Using this parametrization, the circular velocity curve is 
written as: 



d 



dR 



GMR^ 
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(A.4) 



which is independent of the vertical scale length h.. This means 
that any Miyamoto-Nagai model with mass M and horizontal 
scale length Iir will fit the same rotation curve. One may also 
recognize in Eq. |A.4| the familiar velocity curve of the Plummer 
model, the latter being a subclass of the Miyamoto-Nagai model, 
with fl = and b - Iir). 

A corollary of Eq. IA.4I is that the rotation curve of any 
Miyamoto-Nagai model can be reproduced by a linear combi- 
nation of two or more Miyamoto-Nagai models, having different 
scale heights and appropriate masses. For example, the potential 
0;^^(/?,z) defined by 



<1>m(^,z) 



(A.5) 



where / e [0, 1] and /i'. 4^ h,, will have the same velocity curve 



Appendix B: Evolution of model N - 3' 

In this appendix we present the evolution of model - 3', 
where all the gas (visible and dark) is assumed to be dissipa- 
tive, computed with the sticky particle scheme (see Section |23] l. 
The sticky parameters used are similar to the ones used for the 
visible gas in model N = 3 and N = 5. The evolution shows that 
at the beginning (f = - f = 0.8 Gyr) the strongest dissipation 
reinforces the spiral constant in the central part of the galaxy. 
But after 1 Gyr the inner disk, which is dynamically very cold, 
becomes strongly unstable. The spirals break into several dense 
clumps and the disk does not look like a real galaxy. This effect 
is even stronger for an = 5 model, which is not displayed here. 
This test shows that the dark matter in the disk cannot be as dis- 
sipative as the visible disk. Otherwise, the disk is unstable and 
fragments in less than a dynamical time. 
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Fig. B.l. Evolution of model N = 3' where the all the gas (visible 
and dark) is assumed to be dissipative. Each image represents the 
surface density of the visible gas. The box size is 200 x 200 kpc. 
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